## 在这个项目中主要存在4个可以直接优化点，分别是：
# 1、迭代次数 N_sim，如果修改 N_sim 则需要相应的改变图的位置。
# 2、fitness, 表达式或者系数都可以调整。
# 3、突变比例
getwd()
rm(list = ls())
seed <- 3
set.seed(seed)
# install.packages('ggplot2')

library(mgcv) 
library(ggplot2)
library(reshape2) 
library(tidyverse)
source("toolbox.R")


######################################
### settings for simulation number ###
######################################

# parameter
mech <- 5 # choose 1,2,3,4 for other mechanism
N_sim <- 5000  # 迭代次数

# a_list <- c(0,0.1,0.2,0.3,0.35,0.4,0.5,0.6,0.7,0.8,0.9)
# a_list <- c(0.31,0.32,0.33,0.34,0.35,0.36,0.37,0.38,0.39,0.4,0.41,0.42)
# alpha_list <- c(-0.5,-0.3,0,0.3,0.5)
# beta_list <- c(-0.5,-0.3,0,0.3,0.5)

# df <- expand_grid(a_list,alpha_list,beta_list)

# (nrow(df)/5*4+1):(nrow(df)/5*5))

i <- 2

a <- 0.3
b <- 1-a
alpha <- 0
beta <- 0

print(c(a,b,alpha,beta))

# 描述不同fitness定义方式
# for antagonist
F_antagonist <- function(mech){
  E1 <- (1-H_V_A) # n=1 for A_min & P_max
  E2 <- (1-H_V_A)-alpha*(1-H_V_Po)
  E3 <- (1-H_V_A)
  E4 <- (1-H_V_A)+beta*(1-H_V_Po)
  E5 <- (1-H_V_A)+beta*(1-H_V_Po)
  
  E <- list(E1,E2,E3,E4,E5)
  return(E[[mech]])
}

# for pollinators
F_pollinators <- function(mech){
  E1 <- (1-H_V_Po)   # n=1 for A_min & P_max
  E2 <- (1-H_V_Po)-alpha*(1-H_V_A)
  E3 <- (1-H_V_Po)-alpha*(1-H_V_A)
  E4 <- (1-H_V_Po)-alpha*(1-H_V_A)
  E5 <- (1-H_V_Po)+alpha*(1-H_V_A)
  
  E <- list(E1,E2,E3,E4,E5)
  return(E[[mech]])
}

# for plant
F_plant <- function(mech){
  E1 <- (a*H_A_V+b*(1-H_Po_V))
  E2 <- (a*H_A_V+b*(1-H_Po_V))
  E3 <- (a*H_A_V+b*(1-H_Po_V))
  E4 <- (a*H_A_V+b*(1-H_Po_V))
  E5 <- (a*H_A_V+b*(1-H_Po_V))
  
  E <- list(E1,E2,E3,E4,E5)
  return(E[[mech]])
}

###################################################
## load field observation matrix AP_obs, PV_obs ###
###################################################

# # data 应为实际观测值
# AP_obs  <- matrix(data = NA,nrow = 48,ncol = 11) # 这里设置植食昆虫有50种，榕树45种，有机化合物244种
# PoP_obs <- matrix(data = NA,nrow = 18,ncol = 11)
# PV_obs  <- matrix(data = NA,nrow = 11,ncol = 22)

AP_obs  <- read.csv("data/AP_Africa_B.csv", header = TRUE, as.is = TRUE, row.names = 1)
PoP_obs <- read.csv("data/PoP_Africa.csv", header = TRUE, as.is = TRUE, row.names = 1)
PV_obs  <- read.csv("data/PV_Africa.csv", header = TRUE, as.is = TRUE, row.names = 1)


############################
### Preparation:Sample #####
############################

index_choose <- c(35,40,16,40,74,12,14,15,1,19) # 按照原则一抽取的 n 种 VOCs
PV_obs_1 <- PV_obs[,c(index_choose)]
index_all <- 1:ncol(PV_obs)

PV_obs_2 <- PV_obs[,c(index_all[-index_choose])]
PV_obs_3 <- PV_obs_2[,c(colSums(PV_obs[,index_all[-index_choose]])<=i)]

index_1 <- 1:ncol(PV_obs_3)
# index_ <- sample(index_1 , 56)
# print(index_)
# PV_obs_4 <- PV_obs_3[,c(index_ )]

PV_obs <- cbind(PV_obs_1,PV_obs_3)

########################################

nA  <- nrow(AP_obs) # 植食动物的数量 （行计数）
nP  <- ncol(AP_obs) # 植物的数量 （列计数）
nV  <- ncol(PV_obs) #  VOC 的数量
nPo <- nrow(PoP_obs) # 传粉蜂的数量

## 生成模拟的初始化矩阵
# 1、初始矩阵不影响模拟结果
# 2、初始时要保证任意植物至少有与一种植食动物有关；任意VOC都至少由一种植物释放。因此下述的colSum 的所以值都大于 0
AP  <- matrix(rbinom(nA*nP,1,0.5),nA,nP)
PV  <- matrix(rbinom(nP*nV,1,0.5),nP,nV)
PoP <- matrix(rbinom(nPo*nP,1,0.5),nPo,nP)
colSums(AP) #  确保 colSum 的所以值都大于 0
colSums(PV) 
colSums(PoP)
# 每一次循环基因变异的比例，该值仅影响收敛速度，不影响结果。
M1 <- 0.2*nP*nV/2 # PV 矩阵突变的比例
M2 <- 0.2*nA*nP # AP 矩阵突变的比例
M3 <- 0.2*nPo*nP # PoP 矩阵突变的比例

#######################################
### Analyses1: Simulation process #####
#######################################
# create plot
# par(mar=c(4,4,2,9))
png(
  filename = paste0("grid_search_output/date_12_17_seed_",seed,"_",i,".png"), # 文件名称
  width = 650,           # 宽
  height = 480,          # 高
  units = "px",          # 单位
  bg = "white",          # 背景颜色
  res = 72)              # 分辨率

par(mar=c(4,4,2,9))
plot(-1, xlim = c(0,N_sim), ylim = c(0,1), ylab = "Fitness", xlab = "Time")
legend(x=N_sim+N_sim/40, y=0.4, title = "Simulated", legend = c("Plant Fitness","Antagonist Fitness","pollinators Fitness"),
       pch = c(3,4,5), col = c("red","red","red"), xpd=TRUE, bty="n",title.font = 2) 

# legend(x=N_sim+N_sim/40, y=0.8 , title = "Parameter", legend = c(paste0("a = ",a), paste0("b = ",b), paste0("alpha = ",alpha), paste0("beta = ",beta)),
#        xpd=TRUE, bty="n",title.font = 2) 


# 根据实地观察，计算并绘制熵值和fitness
{
  AV_obs  <- as.matrix(AP_obs) %*% as.matrix(PV_obs)
  PoV_obs <- as.matrix(PoP_obs) %*% as.matrix(PV_obs)
  # "H_A_VOC" 函数 i:矩阵 o:互信息/条件熵，toolbox
  # H_V   <- H_A_VOC(PV_obs)[["Hn_V"]]
  # H_P_V <- H_A_VOC(PV_obs)[["Hn_S_V"]]
  # H_V_P <- H_A_VOC(PV_obs)[["Hn_V_S"]]
  H_A_V <- H_A_VOC(AV_obs)[["Hn_S_V"]]
  H_V_A <- H_A_VOC(AV_obs)[["Hn_V_S"]]
  # H_A_P <- H_A_VOC(AP_obs)[["Hn_S_V"]]
  # H_P_A <- H_A_VOC(AP_obs)[["Hn_V_S"]]
  H_Po_V <- H_A_VOC(PoV_obs)[["Hn_S_V"]]
  H_V_Po <- H_A_VOC(PoV_obs)[["Hn_V_S"]]
  # H_Po_P <- H_A_VOC(PoP_obs)[["Hn_S_V"]]
  # H_P_Po <- H_A_VOC(PoP_obs)[["Hn_V_S"]]
  
  # 对抗系统中的 fitness 表达, toolbox
  E_plant       <- F_plant(mech)
  E_antagonist  <- F_antagonist(mech)
  E_pollinators <- F_pollinators(mech)
  
  
  
  # E_obs <- c(0, E_plant, E_antagonist, E_pollinators, H_V,
  #            H_P_V, H_V_P, H_A_V, H_V_A, H_A_P, H_P_A,H_Po_P,H_P_Po,H_Po_V,H_V_Po)
  # names(E_obs) <- c("N","E_plant", "E_antagonist", "E_pollinators", "H_V", "H_P_V", "H_V_P", "H_A_V", "H_V_A", "H_A_P", "H_P_A","H_Po_P","H_P_Po","H_Po_V","H_V_Po")
  #
  ## 野外实际观测值划线，图A
  # 条件熵 for PV, AV and AP
  # abline(h = H_P_V, lty = "dotted", pch = 8)
  # abline(h = H_A_V, lty = "solid", pch = 8)
  # abline(h = H_A_P, lty = "dashed", pch = 8)
  ## 植物和植食动物的fitness，图B
  abline(h = E_antagonist, col = "red")
  abline(h = E_pollinators, col = "blue")
  abline(h = E_plant, col = "green3") # the same as H(A|V) in mech = 1
}


# 模拟过程
A_PV  <- list()  # 每次模拟后存储 PV矩阵 的 list
A_AP  <- list()  # 每次模拟后存储 AP矩阵 的 list
A_PoP <- list()  # 每次模拟后存储 PoP矩阵 的 list
# E = matrix(NA, nrow = N_sim +1, ncol = 15) # E 用来存储每个模拟之后的变量值 N*14
# colnames(E) <- c("N", "E_plant", "E_antagonist", "E_pollinators", "H_V", "H_P_V", "H_V_P", "H_A_V", "H_V_A", "H_A_P", "H_P_A","H_Po_P","H_P_Po","H_Po_V","H_V_Po")
E = matrix(NA, nrow = N_sim +1, ncol = 8) # E 用来存储每个模拟之后的变量值 N*14
colnames(E) <- c("N", "E_plant", "E_antagonist", "E_pollinators",  "H_A_V", "H_V_A","H_Po_V","H_V_Po")


## 基于随机初始化的第一次模拟值
n = 1
{
  AV  <- AP %*% PV
  PoV <- PoP %*% PV
  # "H_A_VOC" 在toolbax.R 中输入矩阵，输出互信息和条件熵
  # H_V   <- H_A_VOC(PV)[["Hn_V"]]
  # H_P_V <- H_A_VOC(PV)[["Hn_S_V"]]
  # H_V_P <- H_A_VOC(PV)[["Hn_V_S"]]
  H_A_V <- H_A_VOC(AV)[["Hn_S_V"]]
  H_V_A <- H_A_VOC(AV)[["Hn_V_S"]]
  # H_A_P <- H_A_VOC(AP)[["Hn_S_V"]]
  # H_P_A <- H_A_VOC(AP)[["Hn_V_S"]]
  H_Po_V <- H_A_VOC(PoV)[["Hn_S_V"]]
  H_V_Po <- H_A_VOC(PoV)[["Hn_V_S"]]
  # H_Po_P <- H_A_VOC(PoP)[["Hn_S_V"]]
  # H_P_Po <- H_A_VOC(PoP)[["Hn_V_S"]]
  
  
  # fitness function
  E_plant       <- F_plant(mech) 
  E_antagonist  <- F_antagonist(mech) 
  E_pollinators <- F_pollinators(mech) 
  
  # 绘图
  # points(n,  H_P_V, col = "red", pch = 24) 
  # points(n,  H_A_V, col = "red", pch=  21)
  # points(n,  H_A_P, col = "red", pch=  22) 
  # 图B
  points(n,  E_pollinators, col = "blue", pch = 5)
  points(n,  E_antagonist, col = "red", pch = 4)
  points(n,  E_plant, col = "green3", pch = 3)
  
  #写入 E, A_AP, A_PV 的第N次
  # E[n,] <- c(n, E_plant, E_antagonist, E_pollinators, H_V,
  #            H_P_V, H_V_P, H_A_V, H_V_A, H_A_P, H_P_A,H_Po_P,H_P_Po,H_Po_V,H_V_Po)
  E[n,] <- c(n, E_plant, E_antagonist, E_pollinators,
             H_A_V, H_V_A, H_Po_V, H_V_Po)
  
  A_AP[[n]] <- AP
  A_PV[[n]] <- PV
  A_PoP[[n]] <- PoP
}


### 植食动物和植物依次使其 fitness 最大化 迭代模拟
for (n in 1:(N_sim/4)) # N_sim : 模拟次数。由于模拟过程包括植食动物、传粉蜂的基因突变和植物两次基因突变四个部分，因此需要 "/4"
{
  print(n)
  ## 传粉蜂：最大化fitness
  for (m in 1:M3) {
    PoP_new <- random.sample_1element1(PoP) #function details see toolbox.r
    PoV <- PoP_new %*% PV
    # H_V   <- H_A_VOC(PV)[["Hn_V"]]
    # H_P_V <- H_A_VOC(PV)[["Hn_S_V"]]
    # H_V_P <- H_A_VOC(PV)[["Hn_V_S"]]
    H_A_V <- H_A_VOC(AV)[["Hn_S_V"]]
    H_V_A <- H_A_VOC(AV)[["Hn_V_S"]]
    # H_A_P <- H_A_VOC(AP)[["Hn_S_V"]]
    # H_P_A <- H_A_VOC(AP)[["Hn_V_S"]]
    H_Po_V <- H_A_VOC(PoV)[["Hn_S_V"]]
    H_V_Po <- H_A_VOC(PoV)[["Hn_V_S"]]
    # H_Po_P <- H_A_VOC(PoP_new)[["Hn_S_V"]]
    # H_P_Po <- H_A_VOC(PoP_new)[["Hn_V_S"]]
    
    
    # fitness function
    E_plant_new       <- F_plant(mech) 
    E_antagonist_new  <- F_antagonist(mech) 
    E_pollinators_new <- F_pollinators(mech) 
    
    # 基于优化机制判断是否保留此次基因突变
    if(E_pollinators_new > E_pollinators)
    {E_antagonis  <- E_antagonist_new
    E_plant       <- E_plant_new
    E_pollinators <- E_pollinators_new
    PoP <- PoP_new 
    }
  }
  
  # 更新矩阵，并计算互信息量和条件熵
  PoV <- PoP %*% PV
  # H_V   <- H_A_VOC(PV)[["Hn_V"]]
  # H_P_V <- H_A_VOC(PV)[["Hn_S_V"]]
  # H_V_P <- H_A_VOC(PV)[["Hn_V_S"]]
  H_A_V <- H_A_VOC(AV)[["Hn_S_V"]]
  H_V_A <- H_A_VOC(AV)[["Hn_V_S"]]
  # H_A_P <- H_A_VOC(AP)[["Hn_S_V"]]
  # H_P_A <- H_A_VOC(AP)[["Hn_V_S"]]
  H_Po_V <- H_A_VOC(PoV)[["Hn_S_V"]]
  H_V_Po <- H_A_VOC(PoV)[["Hn_V_S"]]
  # H_Po_P <- H_A_VOC(PoP)[["Hn_S_V"]]
  # H_P_Po <- H_A_VOC(PoP)[["Hn_V_S"]]
  
  # fitness function
  E_plant       <- F_plant(mech) 
  E_antagonist  <- F_antagonist(mech) 
  E_pollinators <- F_pollinators(mech) 
  
  # 写入模拟完植食动物的 E, A_AP, A_PV
  # E[4*n-2,] <- c(4*n-2, E_plant, E_antagonist, E_pollinators, H_V, 
  #                H_P_V, H_V_P, H_A_V, H_V_A, H_A_P, H_P_A,H_Po_P,H_P_Po,H_Po_V,H_V_Po)
  
  E[4*n-2,] <- c(4*n-2, E_plant, E_antagonist, E_pollinators,
                 H_A_V, H_V_A, H_Po_V, H_V_Po)
  
  A_AP[[4*n-2]] <- AP
  A_PV[[4*n-2]] <- PV
  A_PoP[[4*n-2]] <- PoP
  
  # # 绘图
  # points((4*n),  H_P_V, col = "red", pch = 24) 
  # points((4*n),  H_A_V, col = "red", pch=  21)
  # points((4*n),  H_A_P, col = "red", pch=  22)
  # 图B
  points((4*n),  E_pollinators, col = "blue", pch = 5)
  points((4*n),  E_antagonist, col = "red", pch = 4)
  points((4*n),  E_plant, col = "green3", pch = 3)
  
  ## 被子植物：最大化fitness
  for (m in 1:M1) {
    PV_new <- random.sample_1element1(PV) #function details see toolbox.r
    PoV <- PoP %*% PV_new
    AV  <- AP %*% PV_new
    # H_V   <- H_A_VOC(PV_new)[["Hn_V"]]
    # H_P_V <- H_A_VOC(PV_new)[["Hn_S_V"]]
    # H_V_P <- H_A_VOC(PV_new)[["Hn_V_S"]]
    H_A_V <- H_A_VOC(AV)[["Hn_S_V"]]
    H_V_A <- H_A_VOC(AV)[["Hn_V_S"]]
    # H_A_P <- H_A_VOC(AP)[["Hn_S_V"]]
    # H_P_A <- H_A_VOC(AP)[["Hn_V_S"]]
    H_Po_V <- H_A_VOC(PoV)[["Hn_S_V"]]
    H_V_Po <- H_A_VOC(PoV)[["Hn_V_S"]]
    # H_Po_P <- H_A_VOC(PoP)[["Hn_S_V"]]
    # H_P_Po <- H_A_VOC(PoP)[["Hn_V_S"]]
    
    
    # fitness function
    E_plant_new       <- F_plant(mech) 
    E_antagonist_new  <- F_antagonist(mech) 
    E_pollinators_new <- F_pollinators(mech) 
    
    # 基于优化机制判断是否保留此次基因突变
    if(E_plant_new > E_plant)
    {E_antagonis  <- E_antagonist_new
    E_plant       <- E_plant_new
    E_pollinators <- E_pollinators_new
    PV <- PV_new
    }
    
  }
  
  # 更新矩阵，并计算互信息量和条件熵
  AV  <- AP %*% PV
  PoV <- PoP %*% PV
  # H_V   <- H_A_VOC(PV)[["Hn_V"]]
  # H_P_V <- H_A_VOC(PV)[["Hn_S_V"]]
  # H_V_P <- H_A_VOC(PV)[["Hn_V_S"]]
  H_A_V <- H_A_VOC(AV)[["Hn_S_V"]]
  H_V_A <- H_A_VOC(AV)[["Hn_V_S"]]
  # H_A_P <- H_A_VOC(AP)[["Hn_S_V"]]
  # H_P_A <- H_A_VOC(AP)[["Hn_V_S"]]
  H_Po_V <- H_A_VOC(PoV)[["Hn_S_V"]]
  H_V_Po <- H_A_VOC(PoV)[["Hn_V_S"]]
  # H_Po_P <- H_A_VOC(PoP)[["Hn_S_V"]]
  # H_P_Po <- H_A_VOC(PoP)[["Hn_V_S"]]
  
  # fitness function
  E_plant       <- F_plant(mech) 
  E_antagonist  <- F_antagonist(mech) 
  E_pollinators <- F_pollinators(mech) 
  
  # # 写入模拟完植食动物的 E, A_AP, A_PV
  # E[4*n-1,] <- c(4*n-1, E_plant, E_antagonist, E_pollinators, H_V, 
  #                H_P_V, H_V_P, H_A_V, H_V_A, H_A_P, H_P_A,H_Po_P,H_P_Po,H_Po_V,H_V_Po)
  E[4*n-1,] <- c(4*n-1, E_plant, E_antagonist, E_pollinators,
                 H_A_V, H_V_A, H_Po_V, H_V_Po)
  A_AP[[4*n-1]] <- AP
  A_PV[[4*n-1]] <- PV
  A_PoP[[4*n-1]] <- PoP
  
  # # 绘图
  # points((4*n),  H_P_V, col = "red", pch = 24) 
  # points((4*n),  H_A_V, col = "red", pch=  21)
  # points((4*n),  H_A_P, col = "red", pch=  22)
  # 图B
  points((4*n),  E_pollinators, col = "blue", pch = 5)
  points((4*n),  E_antagonist, col = "red", pch = 4)
  points((4*n),  E_plant, col = "green3", pch = 3)
  
  
  ## 植食动物：最大化fitness
  for (m in 1:M2) {
    AP_new <- random.sample_1element1(AP) #function details see toolbox.r
    AV <- AP_new %*% PV
    # H_V   <- H_A_VOC(PV)[["Hn_V"]]
    # H_P_V <- H_A_VOC(PV)[["Hn_S_V"]]
    # H_V_P <- H_A_VOC(PV)[["Hn_V_S"]]
    H_A_V <- H_A_VOC(AV)[["Hn_S_V"]]
    H_V_A <- H_A_VOC(AV)[["Hn_V_S"]]
    # H_A_P <- H_A_VOC(AP_new)[["Hn_S_V"]]
    # H_P_A <- H_A_VOC(AP_new)[["Hn_V_S"]]
    H_Po_V <- H_A_VOC(PoV)[["Hn_S_V"]]
    H_V_Po <- H_A_VOC(PoV)[["Hn_V_S"]]
    # H_Po_P <- H_A_VOC(PoP)[["Hn_S_V"]]
    # H_P_Po <- H_A_VOC(PoP)[["Hn_V_S"]]
    
    # fitness function
    E_plant_new       <- F_plant(mech) 
    E_antagonist_new  <- F_antagonist(mech) 
    E_pollinators_new <- F_pollinators(mech) 
    
    # 基于优化机制判断是否保留此次基因突变
    if(E_antagonist_new > E_antagonist)
    {E_antagonis  <- E_antagonist_new
    E_plant       <- E_plant_new
    E_pollinators <- E_pollinators_new
    AP <- AP_new 
    }
  }
  
  # 更新矩阵，并计算互信息量和条件熵
  AV  <- AP %*% PV
  # H_V   <- H_A_VOC(PV)[["Hn_V"]]
  # H_P_V <- H_A_VOC(PV)[["Hn_S_V"]]
  # H_V_P <- H_A_VOC(PV)[["Hn_V_S"]]
  H_A_V <- H_A_VOC(AV)[["Hn_S_V"]]
  H_V_A <- H_A_VOC(AV)[["Hn_V_S"]]
  # H_A_P <- H_A_VOC(AP)[["Hn_S_V"]]
  # H_P_A <- H_A_VOC(AP)[["Hn_V_S"]]
  H_Po_V <- H_A_VOC(PoV)[["Hn_S_V"]]
  H_V_Po <- H_A_VOC(PoV)[["Hn_V_S"]]
  # H_Po_P <- H_A_VOC(PoP)[["Hn_S_V"]]
  # H_P_Po <- H_A_VOC(PoP)[["Hn_V_S"]]
  
  # fitness function
  E_plant       <- F_plant(mech) 
  E_antagonist  <- F_antagonist(mech) 
  E_pollinators <- F_pollinators(mech) 
  
  # # 写入模拟完植食动物的 E, A_AP, A_PV
  # E[(4*n),] <- c((4*n), E_plant, E_antagonist, E_pollinators, H_V, 
  #                H_P_V, H_V_P, H_A_V, H_V_A, H_A_P, H_P_A,H_Po_P,H_P_Po,H_Po_V,H_V_Po)
  E[4*n,] <- c(4*n, E_plant, E_antagonist, E_pollinators,
               H_A_V, H_V_A, H_Po_V, H_V_Po)
  A_AP[[(4*n)]] <- AP
  A_PV[[(4*n)]] <- PV
  A_PoP[[4*n]]  <- PoP
  
  # # 绘图
  # points((4*n),  H_P_V, col = "red", pch = 24) 
  # points((4*n),  H_A_V, col = "red", pch=  21)
  # points((4*n),  H_A_P, col = "red", pch=  22)
  # 图B
  points((4*n),  E_pollinators, col = "blue", pch = 5)
  points((4*n),  E_antagonist, col = "red", pch = 4)
  points((4*n),  E_plant, col = "green3", pch = 3)
  
  
  ## 被子植物：最大化fitness
  for (m in 1:M1) {
    PV_new <- random.sample_1element1(PV) #function details see toolbox.r
    PoV <- PoP %*% PV_new
    AV  <- AP %*% PV_new
    # H_V   <- H_A_VOC(PV_new)[["Hn_V"]]
    # H_P_V <- H_A_VOC(PV_new)[["Hn_S_V"]]
    # H_V_P <- H_A_VOC(PV_new)[["Hn_V_S"]]
    H_A_V <- H_A_VOC(AV)[["Hn_S_V"]]
    H_V_A <- H_A_VOC(AV)[["Hn_V_S"]]
    # H_A_P <- H_A_VOC(AP)[["Hn_S_V"]]
    # H_P_A <- H_A_VOC(AP)[["Hn_V_S"]]
    H_Po_V <- H_A_VOC(PoV)[["Hn_S_V"]]
    H_V_Po <- H_A_VOC(PoV)[["Hn_V_S"]]
    # H_Po_P <- H_A_VOC(PoP)[["Hn_S_V"]]
    # H_P_Po <- H_A_VOC(PoP)[["Hn_V_S"]]
    
    
    # fitness function
    E_plant_new       <- F_plant(mech)
    E_antagonist_new  <- F_antagonist(mech)
    E_pollinators_new <- F_pollinators(mech)
    
    # 基于优化机制判断是否保留此次基因突变
    if(E_plant_new > E_plant)
    {E_antagonis  <- E_antagonist_new
    E_plant       <- E_plant_new
    E_pollinators <- E_pollinators_new
    PV <- PV_new
    }
    
  }
  
  # 更新矩阵，并计算互信息量和条件熵
  AV  <- AP %*% PV
  PoV <- PoP %*% PV
  # H_V   <- H_A_VOC(PV)[["Hn_V"]]
  # H_P_V <- H_A_VOC(PV)[["Hn_S_V"]]
  # H_V_P <- H_A_VOC(PV)[["Hn_V_S"]]
  H_A_V <- H_A_VOC(AV)[["Hn_S_V"]]
  H_V_A <- H_A_VOC(AV)[["Hn_V_S"]]
  # H_A_P <- H_A_VOC(AP)[["Hn_S_V"]]
  # H_P_A <- H_A_VOC(AP)[["Hn_V_S"]]
  H_Po_V <- H_A_VOC(PoV)[["Hn_S_V"]]
  H_V_Po <- H_A_VOC(PoV)[["Hn_V_S"]]
  # H_Po_P <- H_A_VOC(PoP)[["Hn_S_V"]]
  # H_P_Po <- H_A_VOC(PoP)[["Hn_V_S"]]
  
  # fitness function
  E_plant       <- F_plant(mech)
  E_antagonist  <- F_antagonist(mech)
  E_pollinators <- F_pollinators(mech)
  
  # # 写入模拟完植食动物的 E, A_AP, A_PV
  # E[4*n+1,] <- c(4*n+1, E_plant, E_antagonist, E_pollinators, H_V,
  #                H_P_V, H_V_P, H_A_V, H_V_A, H_A_P, H_P_A,H_Po_P,H_P_Po,H_Po_V,H_V_Po)
  E[4*n+1,] <- c(4*n+1, E_plant, E_antagonist, E_pollinators,
                 H_A_V, H_V_A, H_Po_V, H_V_Po)
  A_AP[[4*n+1]] <- AP
  A_PV[[4*n+1]] <- PV
  A_PoP[[4*n+1]] <- PoP
  
  # # 绘图
  # points((4*n),  H_P_V, col = "red", pch = 24)
  # points((4*n),  H_A_V, col = "red", pch=  21)
  # points((4*n),  H_A_P, col = "red", pch=  22)
  # 图B
  points((4*n),  E_pollinators, col = "blue", pch = 5)
  points((4*n),  E_antagonis, col = "red", pch = 4)
  points((4*n),  E_plant, col = "green3", pch = 3)
  
}

legend(x=N_sim+N_sim/40, y=0.8 , title = "Fitness", legend = c(paste0("Plant = ",E_plant), paste0("Po = ",E_pollinators), paste0("An = ",E_antagonis)),
       xpd=TRUE, bty="n",title.font = 2) 
dev.off()
# ggsave(
#   filename = paste0("grid_search",i,".png"), # 保存的文件名称。通过后缀来决定生成什么格式的图片
#   width = 7,             # 宽
#   height = 7,            # 高
#   units = "in",          # 单位
#   dpi = 300              # 分辨率DPI
# )




##################################################
##### Analyses2: 随着VOCs数量计算累积相信息  #####
##################################################
# Parameter
N_test <- 31 
Ns = 10000 # set for calculating cumulative mutual information PV
# 生成均匀分布（0~1）的次数，均匀分布用于描述AV这个非01矩阵上的点是否会发生。
# 当AV矩阵上点的概率大于均匀分布时认为这个 AV 这个点描述的事件可以发生，否则则不发生。


########################################
## for herbivore decoding AV matrix
### function -->生成从0到1的均匀分布
M_rdm <- function(r, c) {
  M <- matrix(runif((r*c), 0, 1), nrow = r, ncol = c)
  return(M)
}

M_rdm_AV <- replicate(N_test, M_rdm(nA, nV))


AV <- as.matrix(AP) %*% as.matrix(PV)
Pla <- rowSums(AP)
AV_sc <- AV/Pla

I_AV_result <- NULL
index_list_AV <- c()
for (nr in 1:N_test) {
  print(nr)
  Mr <- M_rdm_AV[[nr]]
  AV <- (AV_sc > Mr)*1 # if the value bigger than random then assign 1, otherwise, assign 0
  I_AV_Ls_r = matrix(NA, nrow = Ns, ncol = nV)
  for(j in 1:nV) {
    print(paste("Library size = ", j, "simulation number =", Ns))
    index_list <- replicate(n= Ns, sample(x = 1:ncol(AV), size = j, replace = FALSE), simplify = FALSE) # 从31个有机化合物中抽取j个
    index_list_AV <- append(index_list_AV, index_list)
    for (i in 1:length(index_list)) {
      P_vi <- Pro_A_Lib(AV[, index_list[[i]]])[[1]]
      IVi  <- Pro_A_Lib(AV[, index_list[[i]]])[[2]]
      I_AV_Ls_r[i,j] <- mut.inf_pr(P_vi, IVi)
    }
  }
  I_AV <- cbind.data.frame(I_AV_Ls_r, "Type" = "AV", "Mech" = "obs", "SimN" = nr) # # change here if it is other mechanism
  I_AV_melt <- melt(I_AV, id.vars = c("Type", "Mech", "SimN"))
  I_AV_result <- rbind.data.frame(I_AV_result, I_AV_melt)
}

  #saveRDS(I_AV_result, paste("I_AV", Mech_vectors[mech],".RDATA",sep = "-"))

I2 <- I_AV_result
colnames(I2)[4]  <- "VOC_number"
colnames(I2)[5]  <- "Mut.Inf"
I2$VOC_number <- as.numeric(I2$VOC_number)

#####
# max for each round of simulation and for each VOC_number
I2_agg1_AV <- I2 %>%
  mutate(index=1:length(index_list_AV)) %>%
  group_by(Mech ,SimN, VOC_number) %>%
  dplyr::summarise(Max.I.Sim = max(Mut.Inf),
                   max.index = index[which.max(Mut.Inf)]) #%>%
index_list_AV[I2_agg1_AV$max.index[2]]
#### due to the non-normal distribution of I_agg_max, take the median value for final analysis
I2_agg_AV <- I2_agg1_AV %>%
  group_by(Mech, VOC_number)  %>%
  dplyr::summarise(Mean.Mut.Inf = median(Max.I.Sim))

# index_list_[13080]

#plot
ggplot(I2_agg_AV, aes(y = Mean.Mut.Inf, x = VOC_number, color = Mech, group = Mech)) +
  geom_line(size = 1, alpha = 0.6) +
  scale_x_continuous(breaks=seq(0, nV, by = 5), limits = c(0, nV)) +
  ggtitle("Herbivore decoding")

## for herbivore decoding PoV matrix
M_rdm_PoV <- replicate(N_test, M_rdm(nPo, nV))


PoV <- as.matrix(PoP) %*% as.matrix(PV)
Pla <- rowSums(PoP)
PoV_sc <- PoV/Pla

I_PoV_result <- NULL
index_list_PoV <- c()
for (nr in 1:N_test) {
  print(nr)
  Mr <- M_rdm_PoV[[nr]]
  PoV <- (PoV_sc > Mr)*1 # if the value bigger than random then assign 1, otherwise, assign 0
  I_PoV_Ls_r = matrix(NA, nrow = Ns, ncol = nV)
  for(j in 1:nV) {
    print(paste("Library size = ", j, "simulation number =", Ns))
    index_list <- replicate(n= Ns, sample(x = 1:ncol(PoV), size = j, replace = FALSE), simplify = FALSE) # 从31个有机化合物中抽取j个
    index_list_PoV <- append(index_list_PoV, index_list)
    for (i in 1:length(index_list)) {
      P_vi <- Pro_A_Lib(PoV[, index_list[[i]]])[[1]]
      IVi  <- Pro_A_Lib(PoV[, index_list[[i]]])[[2]]
      I_PoV_Ls_r[i,j] <- mut.inf_pr(P_vi, IVi)
    }
  }
  I_PoV <- cbind.data.frame(I_PoV_Ls_r, "Type" = "PoV", "Mech" = "obs", "SimN" = nr) # # change here if it is other mechanism
  I_PoV_melt <- melt(I_AV, id.vars = c("Type", "Mech", "SimN"))
  I_PoV_result <- rbind.data.frame(I_PoV_result, I_PoV_melt)
}

#saveRDS(I_AV_result, paste("I_AV", Mech_vectors[mech],".RDATA",sep = "-"))

I2 <- I_PoV_result
colnames(I2)[4]  <- "VOC_number"
colnames(I2)[5]  <- "Mut.Inf"
I2$VOC_number <- as.numeric(I2$VOC_number)

#####
# max for each round of simulation and for each VOC_number
I2_agg1_PoV <- I2 %>%
  mutate(index=1:length(index_list_PoV)) %>%
  group_by(Mech ,SimN, VOC_number) %>%
  dplyr::summarise(Max.I.Sim = max(Mut.Inf),
                   max.index = index[which.max(Mut.Inf)]) #%>%
index_list_PoV[I2_agg1_PoV$max.index[2]]
#### due to the non-normal distribution of I_agg_max, take the median value for final analysis
I2_agg_PoV <- I2_agg1_PoV %>%
  group_by(Mech, VOC_number)  %>%
  dplyr::summarise(Mean.Mut.Inf = median(Max.I.Sim))

# index_list_[13080]

#plot
ggplot(I2_agg_PoV, aes(y = Mean.Mut.Inf, x = VOC_number, color = Mech, group = Mech)) +
  geom_line(size = 1, alpha = 0.6) +
  scale_x_continuous(breaks=seq(0, nV, by = 5), limits = c(0, nV)) +
  ggtitle("Herbivore decoding")

####################
# 互惠对抗1.4版
# 优化代码计算量

